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We analytically calculate the cross sections of double charmonium production in e^e~ — > 
J/tp{'tp{2S))xcj (J=0,l,2) at next-to-leading order (NLO) in Us in nonrelativistic QCD, and confirm 
factorization of these processes. In contrast to XcO production, for which the NLO correction is large 
and positive, the NLO corrections for Xci,2 production can be negative, resulting in decreased K 
factors of 0.91 and 0.78 for J=l and 2 respectively when fi = 2mc. Consequently, the NLO QCD 
corrections markedly enlarge the difference between cross sections of XcO and Xci,2. This may ex- 
plain why e+e" — > J/il){ip{2S))XcO but not e+e~ — ^ J/'tp{'i(){2S))xci,2 is observed experimentally. 
Moreover, for J/tp{'tp{2S))xci,2, the NLO QCD corrections substantially reduce the jj, dependence 
and lead to predictions with small theoretical uncertainties. 

PACS numbers: 13.66.Bc, 12.38.Bx, 14.40.Pq 



I. INTRODUCTION 

The study of heavy quarkonium inclusive produc- 
tion is important to understanding the hadronization of 
heavy quarkonium. In the spirit of nonrelativistic QCD 
(NRQCD) factorization ffl, it is believed that both color- 
singlet and color-octet channels can contribute to heavy 
quarkonium inclusive production. However, due to the 
lack of knowledge of color-octet long-distance matrix el- 
ements, the inclusive production mechanism for heavy 
quarkonium is still not fully understood even though 
complete next-to-leading order (NLO) QCD corrections 
for most of these processes are available |§-0l. On the 
other hand, in heavy quarkonium exclusive production 
all final state particles are targeted, and in some pro- 
cesses the color-octet does not contribute at all, such as 
in double charmonium production in e~^e~ annihilation. 
As a result, the study of double charmonium exclusive 
production, among other kinds of production, may pro- 
vide a good opportunity to learn more about color-singlet 
mechanisms for production and hadronization, without 
theoretical uncertainties from color-octet contributions. 

In 2002, the Belle Collaboration reported a surprisingly 
large double charmonium production cross sectionjpl, 
which was later confirmed by larger data samples |Q, [^ 
and also by BaBar]ll|. The observed double charmonium 
production cross sections for e+e" — ?► J/'ipVciXco) were 
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FIG. 1: Some representative Feynman diagrams for 

e+e" -J> J/-ipXcJ- 



cross section to be larger than that of J/tpVc by a fac- 
tor of 1.8|1^, but experimentally no evidence for dou- 
ble J/ip production was observed. This puzzle was ex- 
plained later by finding a negative NLO QCD correction 
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much larger than the leading order (LO) calculations 
in NRQCD. Aside from attempted explanations by other 
models and methods [p^, in the framework of NRQCD 
these discrepancies were found to be essentially resolved 
by large NLO QCD corrections jlj, [l^l and relativistic 
corrections p^. It is puzzling that, although the cross 
section of J/ipXcO is observed to be very large, J/V'Xci.2 
is not seen. After all, at LO the cross section of J/ipXcO 
is larger than the total cross sections of J/'0Xci,2 by only 
a factor of a little more than 2. Even with the most 
recent dataQi signals of J/ipXci,2 are stiU not clearly 
seen. Recall that there was a similar situation for double 
J/ip production. At LO, NRQCD predicts the J/ij^J/il} 



for double J/jp 
for J/V'Tyc 0|l5 



1 191 and a positive NLO QCD correction 
All these show that the NLO QCD cor- 
rections can be very important for double charmonium 
production. Therefore, it is necessary to examine the 
NLO QCD corrections for J li>Xc\,i as well as J jiixca 
production in e^e~ annihilation and to see whether the 
NLO QCD corrections play important roles in under- 
standing the nonobservation of J/ipXci.2 production at 
B factories. 

The NLO QCD correction for J/ipXcO has been calcu- 
lated in 
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and the difference between XcO and Xci,2 
in the calculation is the summation over polarizations. 
However, the NLO QCD correction for J/ipXci,2 is indeed 
a nontrivial work even though that for J/^XcO has been 



achieved. This is because Xco is a scalar particle, so the 
summation over its spin and angular polarizations can 
be done at the amplitude level, which makes its Lorentz 
structure much simpler than that of Xci 2- However, the 
method used to calculate J/tpXcO in |lq| can hardly be 
adopted here to calculate J/ipXci,2- Furthermore, the 
method in llq] gives only numerical results, while an an- 
alytical expression is really important for analyzing the 
details of the result. Fortunately, using our recently de- 
veloped method, which has been used to do the complete 
NLO corrections for heavy quarkonium hadroproduction 
[^ , analytical expressions of the production cross sections 
of e+e~ — >■ J /il}{i}{2S))xcj at NLO can be conveniently 
achieved. 



II. CALCULATION 

In the following, we briefly describe our calculation. 
We use FeynArts|20| to generate Feynman diagrams and 
Feynman amplitudes. Some representative Feynman dia- 
grams for this process are shown in Fig. |l|. There are gen- 
erally ultraviolet (UV), infrared(IR), and Coulomb singu- 
larities. Conventional dimensional regularization (CDR) 
with D = 4 — 2e is adopted to regularize them. The 
UV- divergences from self-energy and triangle diagrams 
are removed by renormalization. The Coulomb singu- 
lar terms are factored into the J /ij) and XcJ wave func- 
tions. Differing from S-wave charmonium production, 
where IR singularities can be canceled to ensure fac- 
torization, there are generally nonfactorizable IR diver- 
gences in P-wave charmonium exclusive production (such 
as B — ^ XcJ ^10) due to the nonvanishing relative mo- 
mentum between the heavy quark and the antiquark. As 
previously discussed in Q , it is the existence of the as- 
sociated S-wave state J/ip that avoids topologically non- 
factorizable soft interactions and leads to cancellation of 
IR singularities. This argument is generalized in P2| to 
prove the factorization theorem of heavy quarkonium ex- 
clusive production to all orders in as- In the present 
work, we find that the IR singularities are indeed can- 
celed between different diagrams, which further confirms 
factorization of e+e^ — > J /ipXciiJ = 0, 1, 2). 

To perform the calculation, two different methods are 
used, resulting in two completely independent computer 
codes. The results are found to agree with each other 
with high precision. Furthermore, our result for XcO is 
basically consistent with that in [|5| ^. 

One of our calculations is based on the Mathematica 
package FeynCalc[p3|. We separate the soft singular- 
ities in the virtual corrections using the method used 
in p3, and then we treat the singular part analytically 
and the finite part numerically. To calculate the finite 



part, we use the traditional method [g^ to reduce ten- 
sor loop integrations to scalar functions, which are then 
calculated numerically by LoopToolspq. When deriva- 
tives of finite scalar functions are needed, we perform 
them numerically. This method involves large numeri- 
cal cancellations, so we use quadruple precision in the 
calculation to guarantee the result to be reliable. The 
other one of our calculations is based on our self- written 
Mathematica codes. The derivatives are performed be- 
fore tensor reductions and loop integrations, so it avoids 
the derivatives of scalar functions. Then we decompose 
tensor integrals and reduce scalar functions to a funda- 
mental set using integrate-by-part-based recursion rela- 
tions. The fundamental set consists of only one, two, 
and three point scalar functions in this work. We use the 
analytical expressions in |27 for divergent scalar func- 
tions, while finite scalar functions are calculated analyt- 
ically with the methods in (2|]. We use qCDLoop[|^ to 
check our analytical expressions of scalar functions. Fi- 
nal analytical results for the cross sections are given in 
the Appendix. 



III. RESULT AND DISCUSSION 



In the numerical calculation, we set y/s = 10.6 GeV, 
A^ = 338 MeV, and m^ = 1.5 GeV. The values of wave 

MS ' 

functions squared at the origin (or their derivatives) are 
extracted from the leptonic width of J/4> and •4'{2S) and 
the two-photon width of Xc2 129] at NLO in as : 



- Vln-S)— >e+e- 



4a2 / 16a, 



r 



64a^ 



Xc2(n-P)-^77 



i-^-)Kp(o)I 



Setting tts = 0.259 and a ~ j^ we get 

\Ris{0)\^ =0.418GeVm2, 
|^2s(0)P = 0.179GeVm2, 
|i?'ip(0)p = 0.0116GeVTO^. 



(1) 
(2) 



(3) 

(4) 
(5) 



'^ There is an error in the numerical calculation in ]15[ . Our result 
is consistent with it only after correcting the error in ||l5|]. 



Note that the values of wave functions squared at the 
origin here depend on charm quark mass rric. This treat- 
ment can largely cancel the rric dependence in the short- 
distance coefficients and significantly reduce the theoret- 
ical uncertainty due to the choice of rric. 

In Table |l we give the ratio of cross section at NLO 
to that at LO (the K factor) . Different from the XcO pro- 
duction, where the K factor is much larger than 1, the K 
factors for Xci,2 are small, and the NLO corrections are 
small and even negative when fi — 2mc. Therefore, the 
gap of cross sections between XcO and Xci,2 is further en- 
larged at NLO, which gives an explanation for why only 
J/^{'il>{2S))xcO production is observed. In Fig. g, we 
show the cross sections at LO and NLO as functions of 
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FIG. 2: Cross sections of e+e" 
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— > J/ipXcj a-s functions of the renormalization scale ^ at LO and NLO in a^ with 
Vi = 10.6GeV, A = 338MeV, and m^ = 1.5GeV. 




1.45 1.50 1.55 

m,(GeV) 

(a) e+ + e- -> J/ip + XcO 




1.45 1.50 1.55 

/n,(GeV) 



(b)e4 



•^/V' + Xci 



1.6 

1.4 



sS 1-2 









- NLO 1 




' 



1.45 1.50 1.55 

/n,(GeV) 

(c) e+ + e- ^ J/i/j + Xc2 



FIG. 3: Cross sections of e+e — > J/ipXcJ as functions of the charm quark mass tUc at LO and NLO in as with 

Vs = lO.eCeV, A = 338MeV, and ^ = 2toc. 



the renormahzation scale fi. We find that, although /i de- 
pendence for XcO is large at both LO and NLO, QCD cor- 
rections substantially reduce the fi dependence for Xci/2- 
The dependence of charm quark mass is shown in Fig. y, 
where we find the dependence for Xci,2 is also weaker 
than that for XcO- Based on Figs. || and ||, we conclude 
that the NLO predictions for Xci,2 production have small 
theoretical uncertainties; therefore, they can be used to 
precisely test the production mechanism when they can 
be measured in future experiment. 



asi/i) J/^/' + XcO J/^ + Xci J/^ + Xc-2 






2mc 0.259 
Vs/2 0.211 



1.57 
1.79 



0.91 
1.25 



0.78 
1.14 



TABLE I: The K factor of our QCD corrections, with 
Vs = lO.eCeV, A = 338MeV, and m^ = l.SGeV. 

The comparison of experimental data with our theo- 
retical predictions is shown in Table O, where we vary 
rric = 1.5 ± 0.1 GeV to estimate the theoretical uncer- 
tainties. Our results are consistent with all data except 
e+e" — ^ ip{2S)xcO- The Belle Collaboration's observa- 
tion that the cross section of e+e^ — > ■iIj{2S)xcO is big- 



ger than that of e^ 



-^ 



J/ipXcO cannot be explained 



in NRQCD at LO in v, the relative velocity of charm 



quark and anticharm quark, because the only difference 
between ip{2S) and J/ip at LO in v is the wave functions 
squared at the origin, which are well estimated. Perhaps 
the relativistic corrections might give some hint, but the 
large errors in experiment should be reduced before any 
definite conclusion can be drawn. 



IV. SUMMARY 

In NRQCD, using an improved method, we get 
compact analytical expressions of double charmonium 
production cross sections of e+e^ — > J/^{il}{2S))xcj 
(J=0,l,2) at NLO in as and LO in v. Moreover, we 
further confirm factorization of these processes. With 
v^ = 10.6GeV, nic — l.SGeV, and ji = 2mc, we find 
that the cross section for XcO is enhanced by a iiT fac- 
tor of 1.57, while the cross sections of Xci.2 are de- 
creased with K factors of 0.91 and 0.78 respectively. 
The large positive NLO correction to XcO and negative 
NLO corrections to Xci,2 markedly enlarge the differ- 
ence between the cross sections of XcO and Xci,2 and 
provide an explanation for the phenomenon that only 
J/tp{'ip{2S))xcO but not J/tp{ip{2S))xci,2 production is 
observed. Considering the substantially reduced ^ de- 
pendence, and also the small rric dependence, the NLO 
predictions for e~^e~ — > J/ip{'ip{2S))xci.2 are more pre- 
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TABLE II: Comparison of our predicted cross sections with experiments at B factories in units of fb. The errors of 

our theoretical predictions are only from variations of the charm quark mass rric = 1.5 ± O.lGeV. The experimental 

data are cross sections times branchin g fr actions for Xcj decay into more than 2 charged tracks, while the Belle data 

of V-'(2'5') + Xcj in 113 correspond to Xcj decay into at least 1 charged track. 



cise than for many other processes and may be used to 
test the heavy quarkonium production mechanism in the 
future. However, the predicted cross section for 'ip{2S)xcO 
production is much smaller than the data. In order to 
clarify this problem, we suggest that both more careful 
measurement be performed for this process and further 
theoretical investigations be made in the future. 



The LO cross sections for e+ 



J/i' + XcJ are 



a!^° ,= 



hTra^ai 



where A/j° foi' XcO, Xa, and Xc2 are given by 



(A.2) 
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If LO . 


= 32a^ ± 160a^ - 376a^ ± 154a^ ± 45. 
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The NLO cross sections are 



Appendix: Analytical Result 

Here we present the analytical result of our calculation. 
For brevity of expression, we define 
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1, 



l,/ = 8a'* + 4(2a2±l)a5+l, 
1 ^ =8a''-8a^-4a6e±l. 



(A.l) 



where My^o (only real parts are guaranteed to be cor- 
rect) are given by 



^0 = - 24(2560^2 + 78400^° - 14284a* + 8340a*^ - 1943a'' + 69a^ + 30)a^d^li + 96(32a* - 1420a^ 
+ 1710a'' - 321a2 - 45)a'^d^el2 - 24(64ai° + 2904a* + 6360*^ - 784a'' + 420^ + 9)a^(fel3 

- 216(40a'' - 4a2 - 3)a^(fU + 12(256a^2 ^ i9328ai" - 25596a* + 10840a*^ - 1537a'' - 45a2 
+ 21)0^^3^5 + 12(193280^° - 33644a* + 17748a^ - 2867a'' - llla^ + 30)a^d^k + 24(96a^° 

+ 768a* - 812a*^ + 430a'' - 660^ - 9)a^d^eh + 4(l3824ai2 _^ 273152a^° - 152936a* + 20260a^ 

+ 43760" - 1059a2 - 108)aM2e + 2(l92a^° + 11200a* - 13254a^ + 38830" - 810^ - 48)0^^3 ejr^ (A.7) 

- 8(6002688a^^ + 1326848a^" - 2958336a* - 147576a^ + 198772a" + 315240^ + 1377)0^^6 ln2 

- 8(3035136a^2 + 31544320^° - 2263296a* - 276912a^ + 133688a" + 244560^ + 1125)a^6delna 

- 16(80^2 + 656ai" - 2040a* + 4108a*^ - 1796a" + 1260^ + 27)d-''eln(a + b) + 16(ll752a* 
+ 5628a6 - 8364a" + 14300^ + 147)a26c(i3eln(2a + c) + 600(l6a* + 728a^ - 428a" + 380^ 
+ 9)a^bd^e\n{p/mc), 



^NLO ^ ^ 36(896a* + 160a*^ - 1200a" + 1250^ + 60)a^d^eli + 144(848a*^ + 358a" - 4710^ - 8l)a^d^el2 

- 72(l648a* - 5500^ - 436a" + 180^ + 9)02^^6/3 + 18(704a* - 19760^ + 1270a" + 730^ 
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- 48(ll2a^" - 780a* + 1173a*^ - 544a" + 450^ + 27)d^e ln(a + b) + 48(384a* - 5920*^ - 946a" 
+ lOlOa^ + 147)a26cd3eln(2a + c) + 1800(64a^ - 96a" + 260^ + 9)a^bd^e\n{ii/mc), 
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+ 72119360^ + I6O432O0" + 1566600^ + 5625)0^66 In o + 16(17120^^ ^ loiea^" - 3984a* + 469o^ 

+ 1264a" - 387a2 - 135)^^6 ln(a + b) + 16(l808a* + 3192o*^ + 5286o" + 44860^ 

+ 735)a26cd3eln(2o + c) + 600(32a* + 160a^ - 376o" + 1540^ + 45)a^bd^e ln(/x/mc), 



(A.9) 



where h-l^j are combinations of dilogarithms: 



h = Li2 (1 + b/a) , h = Lia (1 - b/a) , h = Li2(e + 2ab) - L^ie - 2ab), 
h - 2Li2(5) - Li2 {dg/f) + Li2 (l//i) - Lia {g/h) - U^ih) + Lis {dh/ f) , 
(a + 6)2] ^. [(a -6)2] ^. [(0 + 6)2] r(a-6)2- 
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